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Abstract This paper presents an adaptive multiple-shooting method to solve 
stochastic multi-point boundary value problems. The heuristic to choose the 
shooting points is based on separating the effects of drift and diffusion terms 
and comparing the corresponding solution components with a pre-specified 
initial approximation. Having obtained the mesh points, we solve the underly- 
ing stochastic differential equation on each shooting interval with a first-order 
strongly-convergent stochastic Runge-Kutta method. We illustrate the effec- 
tiveness of this approach on 1-dimentional and 2-dimentional test problems and 
compare our results with other non-adaptive alternative techniques proposed 
in the literature. 
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1 Introduction 

Numerical methods for solving initial value problems in stochastic differential equations 
(SDE-IVPs) have been extensively researched in the last two decades (see e.g. |17l [20] 
and the references therein). This is not the case for stochastic boundary value problems 
(SDE-BVPs or SBVPs for short), because of complications both in theoretical as well as 
computational aspects. These equations appear naturally in a variety of fields such as 
smoothing [26], maximum a posteriori estimation of trajectories of diffusions [30], wave 
motion in random media [2] , stochastic optimal control |31] , valuation of boundary-linked 
assets [11] and in the study of reciprocal processes [18] . They also arise from the semi- 
discretization in space of stochastic partial differential equations by the method of lines 
approximation [19] . Taking into account the fact that the exact solution of these equations 
are rarely available in analytic form, trying to find efficient approximation schemes for 
the trajectories of the solution process or its moments, seems to be a natural candidate. 
During the last years, several authors have studied with different techniques, the numerical 
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solution of SBVPs of the form: 

/ dX(t) = f(X(t),t)dt + g(X{t),t) odW(t), X(t)eR d , < t < T, 
\a(X)=c, [LL) 

in which / : R d x [0, T] -> R d and g : R d x [0,T] -> M dxd are continuous globally 
Lipschitz functions with polynomial growth, W(t) is a <i-dimensional Wiener process, 
a : C°(R d x [0, T}) — > R d is a continuous operator and c G l rf is a constant vector. The 
existence and uniqueness of the solution process as well as the Markov field property 
of it have been studied by some authors, among them we mention (2H [251 [221 EH ITS'] . 
Due to the anticipative nature of the solution process, the main machinery in the study 
of these equations have turned out to be the Malliavin calculus |21j . 

The majority of research in this field has concentrated around two-point SBVPs (TP- 
SBVPs) corresponding to the choice 

a(X) = h(X(0),X(T))=c, (1.2) 

in which h : R d x R d — >■ R d is a given (possibly nonlinear) function and c is defined as 
before. In this category, we must point out to linear TP-SBVPs in which both the 
drift and diffusion coefficients (/ and g respectively in fll.lj> ) are linear functions of their 
arguments and the function h is of the form 

h(y,z)=H y + H lZ , (1.3) 

in which H$ and Hi are dx d matrices. At the same time, the special class of functional 
boundary conditions of the form 

a{X) = [ dA(t)X(t) = c, (1.4) 
Jo 

have also been of interest, in which A(t) is a d x d matrix valued integrator. The other 
interesting case is the multi-point SBVP (or MP-SBVP for short) having the boundary 
condition 

N s 

a(X) = ^ f A j X(r j ) = c, (1.5) 

in which A\ , A2 ? ■ ■ ■ 5 ^-m 

are constant square matrices of order d and t\, T2, • • • , tn s € [0, T] 
are given switching points with the property Tj < Tj, for i < j. This boundary condition 
could be considered as the result of a quadrature formula applied to approximate the 
general form (|1.4p and will be of special interest in this paper. 

On the numerical side, some efforts have been directed towards devising efficient nu- 
merical schemes for (II. ip among them we mention the following: Allen and Nunn [3] 
propose two methods for linear two-dimensional second order SBVPs, one based on finite 
differences and the other based on simple-shooting. They analyze the convergence proper- 
ties of these methods and report some numerical experiments confirming their theoretical 
results. Arciniega and Allen [5j examine a shooting-type method for systems of linear 
SBVPs of the form (jl.ip . This method could be viewed as a generalization of the comple- 
mentary function approach for deterministic BVPs adopted to solve SBVPs [28] . Arciniega 
[1] extends this work to the nonlinear case and performs some error analysis for this new 
scheme. Ferrante, Kohatsu-Higa and Sanz-Sole |12] use a strong Euler-Maruyama approxi- 
mation to find strong solutions of (jl.ip with linear boundary conditions. They obtain error 
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estimates for this method without accompanying any numerical results to their theoretical 
findings. In a recent paper, Esteban-Bravo and Vidal-Sanz [IT] use the wavelet-collocation 
scheme to find approximations to trajectories of the solution for a general version of (jl.lj) 
with boundary conditions of the form (jl.4p . We must also mention the work of Prigarin 
and Winkler |27j in which they propose a special member of the general Markov chain 
Monte Carlo (MCMC) approach namely the Gibbs sampler to construct realizations of 
the solution process. The convergence is proved for the special case of linear TP-SBVPs 
and some guidelines have been provided to cope with the general nonlinear case and also 
boundary value problems for stochastic partial differential equations. 

Among the above-mentioned schemes, the simple-shooting method which relies on 
transforming the SBVP (jl.ip to an SDE-IVP, has shown to have good accuracy properties, 
but it may give unacceptable approximate solutions on long time intervals. This is specially 
the case when the underlying SDE is unstable i.e. almost all sample paths are rapidly 
growing in absolute value. Our aim here is to circumvent this deficiency by developing 
an adaptive multiple-shooting method to solve (jl.ip based on a detailed analysis of the 
sample paths of the corresponding stochastic equation. The idea is to adaptively subdivide 
the typical interval [Tj,Tj+i] into a grid of shooting points 

in which the tij's and also N(i) will depend on the particular realization (indexed by a;) 
of the underlying Wiener process. In each interval [t^Tj+i], starting from tn = Ti, the 
criterion we choose to obtain ijj+i from ti j is to use an idea adopted from the operator- 
splitting method to investigate the behavior of the two local SDE-IVPs arising from the 
drift and diffusion components of the underlying SDE and controlling upon their growth 
on this subinterval. For this purpose, we employ an initial approximation to the solution 
which (approximately) satisfies the boundary conditions and compare it with the two cor- 
responding SDE-IVP solutions. To obtain the mesh points, we solve the above mentioned 
SDE-IVPs on each shooting interval with a first-order strongly convergent stochastic 
Runge-Kutta method introduced in [13]. We show that this strategy significantly en- 
hances the accuracy and stability properties of the simple-shooting method and at the 
same time reduces the computational cost of the long-time integration problem to a great 
extent. Comparison with other schemes like simple-shooting, finite-differences, wavelet- 
collocation and also the fixed-step multiple shooting method itself, confirms that the pro- 
posed method is a reliable alternative than the widely used non-adaptive approaches in 
the literature. 

The rest of this paper is organized as follows. In section 2, we present the multiple 
shooting framework to solve SBVPs with multi-point boundary conditions. The criterion 
to select the shooting points which forms the foundation of our adaptive strategy will be 
discussed in section 3. The details of optimal parameter tuning for the proposed scheme 
and implementation details will be described in section 4. We conclude the paper by 
commenting on some possible ways to extend this work into more general frameworks. 

2 Multiple Shooting Method for MP-SBVPs 

In this section, we describe the multiple-shooting framework to approximate the sample 
paths of the equation (jl.ip . This can be considered as the extension of methods presented 
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in [5l H] and will serve as the ground base for our adaptive scheme. For this purpose, 
consider the following MP-SBVP in Stratonovich form: 



dX(t) = f(X(t),t)dt + g(X(t),t) o dW(t), X(t) G 



< t < T, 



(2.1) 



Without loss of generality, we assume throughout the paper that t% = and ttv s = T. 
For each realization of the Wiener process, we are interested in finding the corresponding 
realization of the solution process satisfying (|2.ip . Assume that for each i = 1, 2, • • • ,N S — 
1, Ii = [rj,Tj + i] is subdivided into the shooting intervals [tij,tij + i], j = 1, • • • ,N(i) — 1 
with tj 5 i = Tj and i^jvYt) = T{ + \. The adaptive procedure used to obtain them will be 
discussed in section 3 but in the sequel, we assume that they are known. If Xij(t;sij) 
solves the local SDE-IVPs: 



dXij(t;sij) = f(Xij(t;sij),t)dt + g(Xij(t;sij),t) odW(t), t G fcj, tjj+i], 

Xi,j(tij,Sij) — b i,ji 



(2.2) 



for i = 1, 2, ■ ■ ■ , N s — 1 and j = 1, ■ ■ ■ ,N(i) — 1, augmentation of the local solutions 
Xij(t; Si t j) with imposition of continuity condition at the interior shooting points and 
satisfaction of multi-point boundary conditions at switching points will result in a global 
approximation to X(t). For this purpose, we find the unknown initial conditions Sjj's by 
solving the system of D = d x (^Z^x [N(i) — 1]) + 1 nonlinear equations: 

F(s) = 0, (2.3) 

in which 

° — 1*1,1' ' *1,JV(1)-1' *2,1) >*2,jV(2)-l> i b N s -l,l) ) b N s -l,N(N s -l)-l' b N a ,l) fc ^ ) 

is the shooting vector and F(s) is given by: 

s 2 - Xi(si) 



F(s) 



&N a - Xat s _i(sat s _i) 
f(s,c) 



(2.4) 



,N S , 



in which = (s ij2 , s J)3 , • • • , s jiN ^_i, Sj+i,i) T for j = 2, 3, 

^i,i(*i,2;s 3 -,i) 

^7,2(^,3; Sj, 2 ) 

x,. 

^•,AT(j)-2(*i,JV(i)-i; s i,JV(i)-2) 

- x i,iV0')-ife'+ 1 . 1 ' s J.-'V(i)-i) 
for j = 1 , • • • , — 1 and 

#(s,c) = AiSi,! + A 2 S 2 ,1 H h AjV a _iSjv a -l,l + ^Af s SAf s ,l - c. 



(2.5) 



(2.6) 



4 



The solution of system (|2.4p , which provides a global refinement of the solution values at 
the gridpoints, is usually done within the framework of a damped-Newton iteration whose 
k-th. iteration s k is of the form 



s" - A fe [ J DF(s fc )]- 1 F(s fe ). 



(2.7) 



In this relation, Afc £ (0, 1] is the relaxation or damping factor and DF(s k ) is the Jacobian 
matrix of F(s) evaluated at the k-th. iteration. 

It can be shown that 



DF(s) 



-r 2 h 



Ai Ao 



AN a -i A Ns 



in which 



■3,1 



5,2 



r i,AT(j) 



(2.8) 



and the components = D s . h Xj^-\(tjj.', Sj t k-l) for each j and are dxd matrices and 

Idxd 



Id> 



L dxd 



(2.9) 



is an N(j) x N(j) identity matrix. It is obvious that the exact computation of Tj^ 
requires the analytic solution of the local SDE-IVPs ()2.2|) . It is worth pointing out here 
that although it is possible to approximate I^'s by linearization of the corresponding 
local SDEs and integrating them up to tij, we will adopt an alternative strategy by 
approximating the derivative terms by finite differences. A strategy for choosing the Afc's 
has also been developed and thoroughly tested in [9] which will be pursued here. 



3 Adaptive Sequential Selection of Shooting Points 

Multiple-shooting method as a natural generalization of the simple-shooting idea, signif- 
icantly enhances the stability properties of its ancestor and behaves much better than it 
in terms of accuracy and rate of convergence and so has been a preferred choice to solve 
deterministic boundary value problems [HJ \TE\ 129] . The main drawback of this method 
could be attributed to its computational cost which is directly proportional to the number 
of shooting points in the integration interval. To reduce these costs, some authors have 
proposed to devise a control mechanism on the number and location of shooting points in 
such a way that the stability and accuracy properties of the method are preserved. This 
strategy has the additional advantage of resolving the special features of the solution in 
the integration interval: "... a multiple-shooting approach should permit step sizes to be 
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chosen sequentially, fine in the boundary layers, and coarse in the smooth regions" |10j . 
We extend this argument to the case of non-smooth solutions - the feature which is intrin- 
sic for SDEs - and show that the adaptive selection of shooting points based on the driving 
force for this non-smooth behavior, i.e. the underlying Wiener process and also comparing 
the solution with an initial approximate solution, will have an overall performance much 
better than the corresponding fixed step-size counterpart. 

To start the adaptive procedure, we first find a simple piecewise linear approximation 
to the solution, 0(t), which approximately satisfies the multi-point boundary conditions 
(|1.5p . To find this approximation, we discretize the SDE-IVP in each interval with the 
Euler-Maruyama method and then solve the following system of nonlinear equations for 
" 's: 



9_ T2 = 9 T1 + (r 2 - ri)/(5 n ,ri) + (W(r 2 ) - W(n)) 5 (0 n , r x ) 
6t3 = T2 + ( T3 _ r 2 )f(9 T2 ,r 2 ) + (W(t 3 ) - W(T 2 ))g(9 T2 ,r 2 ) 



'r Ns "t Ns 



-1 + ( T N S ~ TN.-l)f(B T jt a _ 1 ,TN,-l) + (W(t Ns ) - W(TN,-l))g(Ont,- J .,TN.-i), 



(3-1) 

The continuous piecewise linear approximation could then be obtained by linear interpo- 
lation: 

o(t)= — e n+1 + Ti+1 ~ t 9 Ti , te[T h T i+1 ], * = i, - - - ,n s . (3.2) 
n+i — n Tj+i — Ti 

Consider now the interval [rj,Tj_|_i] and put i^i := Tj. Starting from tij and to obtain the 
next shooting point in this interval, we integrate the following two local SDE-IVP's: 

f dX(t) = f(X(t),t)dt, te[t itj ,T i+1 ], (33) 
\ x(uj) = o(ti,j), 1 ' ; 

f dX{t) = g (x(t),t)odW(t), te[t itj , n+1 ], 

\ X(t hj ) = 9(tij), 

by deterministic and stochastic components of an SRK method, described in the next sec- 
tion. We will terminate the integration when we reach the first point in our discretization 
satisfying: 

t id <s<n + x, \\X(s)\\>Li(8) or \\X(s)\\>L 2 {s), (3.5) 

in which L\(s) and L 2 {s) will be specified in the sequel. We then put tij+i := s and restart 
the integration of both (|3,3p and f|3.4[) from tij+i using X (tjj+i) = X(tij + i) = 8(tij + i) 
as the initial guess. This procedure will be continued up until the point f$ Mt) = r «+i i s 
reached and then will be continued from the next shooting interval to finally arrive at T. 

The first untold story in our description of the algorithm is the selection of the "stop- 
loss functions" Li(s) and L 2 {s) which control upon the location of our shooting points. 
The most intuitionistic proposal could be 

Li( a ) = a||0( a )||, L 2 ( S ) = /3||0( S )|| (3.6) 



for some positive constants a and j3 (see e.g. [29J Section 7.3.6 for a similar idea in the case 
of deterministic BVPs). We could choose the a and j3 coefficients in (|3,6p time-dependent 
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and find an empirical optimal relation for them, but our numerical experience shows 
that the gain in efficiency is not substantial. We have also tested other stopping criteria 
based only on the size of the increments of the Wiener process which has resulted in the 
selection of more shooting points but has not improved the accuracy in a comprehendible 
way. Another proposal is to find the first point s which simultaneously maximizes the 
following quantities: 



\X(s)\\ ^ \ „f\\X(s) 
> a 



\\9(s)\\ ~ r \\\9{s)\\ 



>P) 



for given positive a and (3. The idea has led us to solve simple constrained stochastic 
programming problems in each step (with exact solutions for linear SBVPs) that needs 
further investigation and will be pursued in a forthcoming paper. 

It is evident from the form of our adaptation criteria that in the case of weak driving 
noise process, we are controlling upon the size of the solution process and look at the first 
time at which the norm of the solution starts to deviate from the initial piecewise linear 
approximation. On the other hand, when the increments of the Brownian noise become 
large in some portions of the solution domain, we must finish the integration and select 
the current point as a suitable shooting point. In both of these scenarios, we must come 
back to the initial approximation 9(t) and continue the integration from the initial value 
Q(U,j+i)- 

Having described the way in which we choose our shooting points for each realization, 
we now need to tell the other story about the time marching procedure to solve our SDE- 
IVPs resulting from the multiple shooting method, which will be discussed in the next 
section. 



3.1 Stochastic Runge-Kutta Family 

Among the many possible choices of methods to integrate the ODE-IVP and SDE-IVP 
problems in (13.3j) and (|3.4p (assuming w.l.o.g. that both equations are autonomous), we 
choose to work with a special member from the general class of stochastic Runge-Kutta 
(SRK) methods of the form 

Vi = x n + hY,j=iaijf(Vj) + JiT,j=ibij9(Vj), i = l,...,s 

(3-7) 

x n+1 = x n + /iEi=i a jf(vj) + Ji Ei=i ijdivj), 

in which X n and X n+ \ are approximations to X(t n ) and X(t n+ \) respectively, h = t n +i—t n 
and Jx = odW(s) = W(t n+1 ) - W(t n ). Here, A = (c%-)^- = i and B = (^)™ J=1 are 
s x s matrices with real elements and a T = (ai, . . . ,a s ) and ^ T = (71, . . . ,j s ) are row 
vectors in R s . A typical member of this family could be represented by the Butcher tableau 





A 


B 




T 

a 


T 

7 



and according to the theory presented in [7], the highest possible order of strong (and 
also weak) convergence among all consistent choices for A,B,a and 7 is one (see [T7] for 
notions of strong and weak convergence in the SDE literature). In this work, we use a 
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three stage SRK method (dubbed R3) as the underlying numerical integrator which has 
the tableau 

























1 

2 








1 

2 













3 
4 








3 
4 







•2 


3 


4 


•2 


3 


4 




9 


9 


9 


9 


9 


9 



and its deterministic and stochastic components are themselves valid numerical integration 
schemes [13] . More specifically, we integrate the IVP (|3.3|) with a method of the form 



V2 

X n+1 = X n + |(2/(^) + 3/(ry 2 ) + 4/(r7 3 )), 



(3.8) 



and integrate the SDE (I3.4p with another method having the form 



m 

x n ^ 



Xr 



x n + %g(vi), 

X n + ^g(rj 2 ), 
X n + f(2g(m) 



(3.9) 



+ 3g(rj 2 )+4g(fj 3 )). 



It is interesting to note here that the first scheme has third-order of convergence for a 
deterministic IVP and this will result in higher precision when we are faced with an SDE- 
BVP having a weak driving noise. On the other hand and for the second scheme, we have 
first-order of strong convergence for drift-free SDEs and when the drift is going to diminish 
in some portions of the problem domain, we have an exact-enough method to trace the 
non-smooth path of the corresponding realization. 



4 Numerical Experiments 

In this section, we report on the numerical results obtained using the adaptive multiple- 
shooting method proposed in this paper. We compare its performance with that of its 
peers, namely a method based on wavelet-collocation introduced in [IT], a finite-difference 
scheme first analyzed in [3] and adopted here to solve multi-point SBVPs (see the Appendix 
for details of its derivation) and a simple-shooting method when it applies. 

We have selected three test problems from the literature each exemplifying differ- 
ent characteristics of the solution process. The first problem is a 1-dimensional SBVP 
with a functional boundary condition and additive noise but the other two are linear 2- 
dimensional TP-SBVPs, the first with additive and the second with multiplicative noise. 
All of the algorithms are implemented in the MATLAB problem-solving environment and 
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executed on a core i5 processor, 2.4GHz, 4GB RAM computer. 



Test Problem 1: In this numerical experiment, we try to solve the following 1-dimensional 
SBVP with functional boundary condition 

f dX(t) = lo dWCt), < t < 1, 

{ i ~ 

and having the exact solution [TT] 

X(t) = - [ W{s)ds + W t . (4.2) 
Jo 

The integral boundary condition in (|4.ip should be discretized (e.g. by the trapezoidal 
method) into a multi-point boundary condition of the form 

A Ns_1 A 

-^x( n ) + Y, ArXfo) + -Y x i7N a ) = o, 

i=2 

in which Tj = (j — l)Ar for j = 1, • • ■ , N s is the j-th. switching point and Ar = N 1 _ 1 . We 
now place iV m equally-spaced points on the interval = [Tj,T i+1 ] which act as the base 
mesh to integrate the resulting local SDE-IVPs and the global SDE problem. For each 
realization of the Wiener process (constructed on the base mesh), we solve the system 
of equations (|3.ip for 9 Tj , j = 1,2, •• • , N s and interpolate them by (|3.2p to arrive at a 
globally-defined piecewise linear initial approximation to the solution on the whole unit 
interval. 



To find the location of shooting points on /j, we start to synchronously integrate (|3.3|) 
and (|3.4p on the base mesh with the schemes described in Section 3.1 to arrive at the 
first point satisfying (|3.5p with a = and /? = 2.5, from where we turn back to 9{t) and 
continue the process to reach Tj+i. Similar procedure must be repeated for other intervals 
to find the set of all optimal shooting points on [0, 1]. 

We are now ready to form and solve the nonlinear system (|2.3p by a damped-Newton 
iteration (adopted from [9]) to obtain the optimal starting values at the shooting points 
(sj, i = 1, 2, • • • , N s in (I2.4p ) and finally solve the original SDE with these initial values 
by the underlying (full) R3 scheme. 

The accuracy of different schemes is measured via the measure which is an average 
of the form 



k=l 

over a fixed number M of realizations from the maximum grid-wise error 

E(u k ) := max \X(ti,u k ) - X(ti,u k )\ 

1=1,2, ■■■ ,N 

approximating the expected supremum norm 

E(\\X(t)-X(t)\\ c 
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in which X(t) is our approximation and X{t) is the exact solution. We have used the 
quadl function in MATLAB to approximate the integral term in the exact solution (|4.2p 
which uses the adaptive Lobatto quadrature method. 

The results of our computations are depicted in Tables 1 and 2. To compare the 
accuracy over a single realization (M = 1), we have provided Table 1 with columns 
reporting the global error (Eoo) for two different methods and a range of grid spacings 
in the problem domain. The meaning of N in the wavelet-based method is the number 
of collocation points and in the adaptive multiple-shooting method (or adaptive MSM for 
short) is the size of base grid used in the integration process. It must be noted that the 
number of switching points (N s ) we have used in the i-th row is chosen to be 2 l and the 
number of mid-points (N m ) is set accordingly. The superior accuracy of the proposed 
method (granting one order of magnitude more precision in the results) is obvious from 
this table. We have observed similar patterns of error behavior over many realizations 
[M » 1) for the adaptive multiple-shooting but due to the unavailability of the data for 
the other scheme, we have not included them in the Table HI 

Table 1: Comparing the accuracy of the wavelet-based and adaptive multiple-shooting 
methods. 



N 


Wavelet-Collocation 


Adaptive MSM 


2 2 


0.2058 


0.0266 


2 1 


0.0997 


0.0036 


2 6 


0.0075 


0.0007 



We also have compared the proposed method with a finite-difference scheme in Table 
2. Here the errors are reported over M = 500 realizations (in both methods) and the finite 
difference equations are set up on all of the N grid-points of the base grid in our adaptive 
scheme. The column with the heading N a indicates the average number of shooting points 
selected by the algorithm. Again we observe a higher accuracy for the adaptive method 
and a rapid rate of convergence to the exact solution. 

Table 2: Comparing the accuracy of the finite-difference and adaptive multiple-shooting 
methods. 



N 


N s 


N m 


N a 


FD Method 


Adaptive MSM 


2 b 


7 


4 


11 


0.4819 


0.0377 


2 6 


10 


6 


16 


0.2728 


0.0264 


2 7 


15 


8 


23 


0.2477 


0.0164 


2 8 


22 


12 


36 


0.2409 


0.0111 


2 9 


32 


16 


52 


0.2287 


0.0074 



In order to investigate the rate of decay of the error (in the strong sense) for the adap- 
tive multiple-shooting method, we have plotted Figure 2] which shows, in a logarithmic 
scale, the behavior of the global error in terms of increasing the number of switching points. 
One can observe that the rate of convergence is linear in At and the line of linear regression 
applied to the data has a slope of q = 1.0126 with a residual r = 0.0908. This is a priori 
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- Adaptive MSM Global Error 

- Reference Slop of 1 



Figure 1: Asterisks: strong error measure for the adaptive multiple-shooting method 
applied to test problem (1). Dashed line: reference slope of 1. 

anticipated as we have used a method of strong order of convergence one in the integration 
procedure and a super-linear convergent method in solving the set of nonlinear equations. 



Test Problem 2: Here we solve the following 2-dimensional TP-SBVP 



in which 



and 



dX(t) = (AX(t) + a)dt + (BX(t) +b)o dW u 
H o X(0) + ffiJf (1) = c. 



A 



<t < 1, 



" 


1 " 







, B = 


' 


" 













, a = 










, b = 




c i 






. C 2 _ 



Hn 



' 1 


" 


, H x = 


' 


" 




" " 








1 













Writing X(t) = [X\(t), X2(t)] T , it could be shown that Xi(t) solves a second-order SDE 
having the exact solution 

X 1 (t) = c 1 t -^ + c 2 (t-l)^sdW(s) + tf t 1 (s-l)dW(s), (4.3) 

and X 2 (t) = (see [3] for more details). 

To obtain the initial trajectory 9(t) = [6i(t), 02{t)] T for this test problem and supposing 
that a realization of W(l) is simulated, we first solve the linear system 



H + H 1 + H X A + W(l)H 1 B)0 n =c-H ia - W(l)b 



for 9 T1 and then solve another linear system 



H\6 T2 — c — H[ 



for T2 . Now we use linear interpolation to obtain 9{t) over the whole unit interval. In 
checking the relation (|3.5p . we have used the /oo-norm on both sides with a = 2 and 



11 



(3 = 1.5. We also compute the integral terms in the exact solution (|4.3p by a highly- 
accurate trapezoidal scheme. 

We use the finite-difference and also the simple-shooting methods as two competing 
approaches to solve this same problem. Table 3 summarizes our computational results for 
the case c\ = 1 and c% = 1 averaged over M = 1000 realizations. We can observe that 

Table 3: Comparing the accuracy of the finite-difference, simple-shooting and adaptive 
multiple-shooting methods. 



N 


N a 


FD Method 


SS Method 


Adaptive MSM 


2 5 


9 


0.0121 


0.0123 


1.14e-16 


2 6 


11 


0.0064 


0.0065 


1.70e-16 


2 7 


13 


0.0032 


0.0032 


2.59e-16 


2 8 


15 


0.0016 


0.0016 


4.62e-16 


2 9 


17 


0.0008 


0.0008 


9.62e-16 



while both finite-difference and simple-shooting methods converge uniformly to each other 
(in terms of accuracy and order of convergence), the adaptive multiple-shooting beats 
them and gives very accurate results. We also observe a steady growth in the errors as 
we increase N which could be attributed to the accumulation of round-off errors in the 
solution process. 

In order to show the efficiency of the adaptive method in the weak sense and using the 
fact that we can compute the expectation of the exact solution and its non-central second 
moment by the following formulas 

nxi(t)) = ci^y^, mm = (H+^d^^, (4. 4 ) 

we have approximated these expected values on a range of points in the solution domain 
by averaging over M = 10000 realizations of the solution process (computed pointwise) 
and have compared the results with that of other schemes listed in Table 4. While all 
methods have a comparable accuracy, the performance of the adaptive method is actually 
slightly better at all points in the range [0, 1]. 



Table 4: Comparing the accuracy in the weak-sense of the Heun simple-shooting, finite- 
difference and adaptive multiple-shooting methods. 





Heun Simple- 


-Shooting [3] 


FD Method [3] 


R3 Adaptive MSM 


Exact 


t 


E{X{t)) 


E(X*(t)) 


E(X(t)) 


E{X\t)) 


E(X(t)) 


E(X 2 (t)) 


E(X(t)) 


E(X 2 (t)) 


0.0 


-0.0000 


0.00000 


-0.0000 


0.00000 


-0.0000 


0.0000 


-0.0000 


0.0000 


0.2 


-0.0800 


0.01499 


-0.0805 


0.01497 


-0.0800 


0.0151 


-0.0800 


0.0149 


0.4 


-0.1194 


0.03338 


-0.1201 


0.03357 


-0.1199 


0.0338 


-0.1200 


0.0336 


0.6 


-0.1192 


0.03346 


-0.1193 


0.03328 


-0.1202 


0.0338 


-0.1200 


0.0336 


0.8 


-0.0793 


0.01486 


-0.0791 


0.01472 


-0.0800 


0.0150 


-0.0800 


0.0149 


0.1 


-0.0000 


0.00000 


-0.0000 


0.00000 


-0.0000 


0.0000 


-0.0000 


0.0000 



Test Problem 3: As the last example, we solve the 2-dimensional SDE-BVP system 
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(adopted from [23]) of the form: 

( dX(t) = BxXit) o dWx(t) + B 2 X{t) o dW 2 (t), < t < 1, 
1 F X(0) + fTiX(l) = c, 



in which 



and 



Si 



1 1 




, B 2 




1 



Hn 



' 1 


1 " 


, H x = 


' 


" 




' 1 " 











1 _ 


, c = 


1 





This equation has an exact solution of the form 



where 



X(t) 



e w 1 (t)( l _ e -W 2 (i) +a o e -W 2 (i) 
P W 2 (t) _ p Ws(l) 



(4.5) 



(4.6) 



Similar to test problem (2), we could obtain the initial trajectory 9(t) = [Oi(t), 02{t)] T by 
first simulating a realization from W(l) = [Wi(l), W^l)] 7 * an d then solving the two linear 
systems 



(H + Ht(I + BxW^l) + B 2 W a (l)))0 n 

H\9 T2 



c-H, 



(4.7) 
(4.8) 



for 9 T1 and T2 respectively. Now 8(t) is computed by linear interpolation and the inte- 
gration is started to obtain the location of shooting points in the base grid. We use the 
absolute values of the second components of X, X and in (|3.5h with a = 1.5 and /3 = 2 
and approximate the integral terms in (|4.5p and (|4.6p by a sufficiently accurate trapezoidal 
scheme. 

The results of our computations for this test problem are reported in Table 5. For 
comparison purposes, we have also included the results of applying fixed-step multiple- 
shooting method in this table. To be fair in the competition, we have selected the number 
of shooting points in the fixed-step multiple-shooting equal to the average number of 
adaptive shooting pointes (N a ) selected by the adaptive algorithm. 

In order to examine the strong order of convergence of the adaptive scheme, we have 
prepared Figure H] which shows clearly (and in a logarithmic scale) that this order is one. 
The result of linear regression applied to the data used in the figure gives us a slope of 
q = 1.0515 with residual r = 0.1190 which is acceptable. 



5 Concluding Remarks 

The numerical solution of boundary value problems in stochastic differential equations is a 
highly unexplored territory of the SDE world requiring the special attention of the experts 
in the field to devise methods of high accuracy and efficiency with low computational 
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Table 5: Comparing the accuracy of the fixed-step multiple-shooting and adaptive 
multiple-shooting methods. 



N 


iV a 


Fixed MSM Method 


Adaptive MSM R3 


2 b 


3 


0.0137 


0.0093 


2 6 


4 


0.0049 


0.0041 


2 7 


4 


0.0025 


0.0021 


2 8 


5 


0.0012 


0.0009 


2 9 


5 


0.0006 


0.0005 




Figure 2: Asterisks: strong error measure for the adaptive multiple-shooting method 
applied to test problem (3). Dashed line: reference slope of 1. 

demand and complexity. We have proposed in this paper, an adaptive multiple-shooting 
method for general multi-point SBVPs based on a stochastic Runge-Kutta integrator. 
Although the adaptation criteria is simple and easily implementable in the method, it 
gives acceptable results in comparison with some other non-adaptive alternatives proposed 
in the literature. The next step in our research (as explained briefly in Section 3) is to 
make use of more elaborate stopping criteria in the selection of shooting points and its 
theortical analysis. We could also incorporate the idea of adptive time-stepping in the 
integration process itself which we anticipate to improve the accuracy further but needs 
a theoretical foundation to prove the stability of the overall scheme in a unified manner. 
Finally, we must also mention the need for introduction of nonlinear test instances into the 
field which is of great importance for testing and benchmarking purposes in the algorithmic 
developements expected to be seen in the near future. 

A Finite-Difference Method for Multi-Point SBVPs 

In this appendix, we present a finite-difference scheme for multi-point SBVPs of the form 

L[X](t) = dW{t) (A.l) 
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in which the operator L is defined by 



L = D n + a n ^(t)D n + ■■■ + ai (t)D l + a (t), D :-- 



where the coefficients aj(i)'s, i = 0, 1, ■ ■ ■ , n — 1 are continuous functions defined on [0, 1]. 
We also append (lA.lj) with boundary conditions of the form 



A'., 



3=1 



1 < i < n. 



(A.2) 



defined on some switching points < Tj < 1, j = 1, 2, • • • , iV s (see [2] for a detailed study 
of some features of the solution to these problems). 

Similar to ordinary differential equations, the SBVP ()A,ip - (|A.2p can be turned into a 
first order system 

dY(t) + A(t)Y(t) = dW(t), (A.3) 



constrained to satisfy 



^ OLijY n (Tj) = a, l<i < 

3=1 



n. 



(A.4) 



in which Y(t) = (Yi (*),••• , Y„(t)), = D^Xit) for 1 < i < n, W(t) = (W(t),0,--- ,0) 
and 



A(t) 



a„_i(t) a n _ 2 (i) 
-1 
-1 











ai{t) a (t) 





1 







To solve (|A.3|) - (|A.4|) by the finite-difference method, we first construct a base mesh in- 
cluding the switching points of the form 

= h < t 2 < ■■ ■ < tjf-i < t N = 1 

and use an explicit one-step difference scheme on this mesh to arrive at 

Y f - : - Y j + A j Y j = AW 3 , 

in which Y J = [Y/, , . . . , Yi] T is an approximation to Y(tj), A> := A(tj) and AW J := 
W(ij + i) — W(tj). Simplifying the above relation we obtain 

Y j+1 + (A j - I)Y j = AW j 

and arranging them in a sequential manner (into a linear structure) we reach to 

AY = w (A.5) 



in which 



Aj = Aj - I, A 



Ai I 

A 2 / 

Atv-i I 

$1 $ 2 ••• $N-1 ®N j {Nxn)x(Nxn) 
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and 





" 


. 


. 


ay 




ay 







. 


. 


a 2j 




a 2 j 




_ 


. 


. 


Oi n j _ 







x [0,0,..., 0,1]! 



Also the vector Y has the form 

Y=[Y 1 ,Y 2 ,...,Y JV ] T 

and 

w = [AW 1 , AW 2 , . . . , AW^f. 

By solving ()A.5h for each realization of the Wiener process, one obtains the corresponding 
realization for the solution process on the base mesh which is what we have reported for 
test problem (1). 
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